################################################################################
#  South-South migration: the impact of the Venezuelan diaspora on Colombian   #
#                     natives' wages and hours worked                          #
#                          Applied Economics Letters                           #
#                    Mora, Cuadros-Menaca, and Sayago (2022)                   #
################################################################################

rm(list=ls()) #start clean

#Packages

dev.off(dev.list()["RStudioGD"])
ipak <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE, repos = "http://cran.us.r-project.org")
  sapply(pkg, require, character.only = TRUE)
}
ipak(c("rstudioapi", "readxl", "tidyverse", "data.table", "janitor","lfe", "sf", "splm", "spatialreg",
       "spsur","plm","openxlsx","RColorBrewer", "dplyr","maptools", "stargazer","sf","spdep", "openxlsx",
       "ggplot2","classInt","cowplot", "googleway", "ggplot2", "ggrepel", "readstata13","gghighlight","ggpubr",
       "ggspatial", "tidytext","stringr","tm",'tidyr',"wordcloud2","textdata","syuzhet","tidyselect","reshape2",
       "data.table","wordcloud", "sf", "rnaturalearth", "rnaturalearthdata" ,"rtweet","academictwitteR",
       "lubridate","readstata13","foreign"))

setwd('Your Directory')

#Shape Files

shape=st_read("Municipios.shp")
shape2=read_sf("Municipios_Venezuela.shp")

#Data

raw1=read.xlsx("tablamapa.xlsx")
raw1$Big[is.na(raw1$Big)]=0
shapemod=merge(shape,raw1,by="ID_ESPACIA")

raw2=read.xlsx("cod_mun.xlsx")
raw2$ID=as.integer(raw2$ID2)
shapemod2=merge(shape2,raw2,by="ID")

#Map

world <- ne_countries(scale = "medium", returnclass = "sf")
class(world)

names(world)
head(world)
table(world$iso_a3)
world=subset(world,continent=="South America")
world2=subset(world,iso_a3=="COL"|iso_a3=="VEN")

centroids.df <- st_coordinates(world2)
world_points<- st_centroid(world2)
world_points <- cbind(world2, st_coordinates(st_centroid(world2$geometry)))
head(world_points)

centroids.df <- st_coordinates(shapemod)
mun_points<- st_centroid(shapemod)
mun_points <- cbind(shapemod, st_coordinates(st_centroid(shapemod$geometry)))

names(shapemod)
shapemod$Treatment2=shapemod$Treatment
shapemod$Treated=factor(shapemod$Treatment)
mun_points$Treatment2=mun_points$Treatment
mun_points$Treated=factor(mun_points$Treatment)

names(mun_points)

breaks_qt <- classIntervals(c(min(mun_points$Temperatura) - .11, mun_points$Temperatura+1), n = 3, style = "jenks",na.omit=T)
breaks_qt

mun_points <- mutate(mun_points, Temperature = cut(Temperatura, breaks_qt$brks)) 

names(shapemod2)
summary(shapemod2$Temperatura.promedio)
shapemod2 <- mutate(shapemod2, Temperature2 = cut(Temperatura.promedio, breaks_qt$brks)) 

 ggplot(data=mun_points) +
   geom_sf(data = world_points,color = "black", fill = "white")+
   geom_text(data= world_points,aes(x=X, y=Y, label=admin),
            color = "black", fontface = "bold",size=2, check_overlap =T) +
   geom_point(data=mun_points,aes(x=X, y=Y,colour=Treated),size=2)  +
   geom_text(data= mun_points,aes(x=X, y=Y, label=NOM_MUNICI.y),
               color = "darkblue", fontface = "bold",size=2, check_overlap = T)+
   theme_minimal()
 ggsave("map_Treated.png",  dpi = "retina")    

ggplot(data=mun_points) +
  geom_sf(data = world_points,color = "black", fill = "white")+
  #  geom_text(data= world_points,aes(x=X, y=Y, label=admin),
  #            color = "black", fontface = "bold",size=2, check_overlap =T) +
  geom_sf(data=mun_points,aes(fill=Temperature))+
  geom_sf(data=shapemod2,aes(fill=Temperature2)) + scale_color_brewer(palette = "Reds") +
  geom_text(data= mun_points,aes(x=X, y=Y, label=NOM_MUNICI.y), color = "darkblue", fontface = "bold",size=2, check_overlap = T)+
  theme_minimal()
ggsave("map_Temperature3.png",  dpi = "retina")    





